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SECTION  I 


GENERAL  DISCUSSION 

The  SAI  program  in  Nuclear  Airblast  Effects  under 
contract  number  N00014-81-C-2267  has  addressed  the  improve¬ 
ment  of  the  numerical  algorithms  and  the  simulation  of  high- 
explosive  and  nuclear  detonations.  Significant  improvement 
in  the  basic  numerical  algorithm,  Flux  Corrected  Transport, 

FCT ,  has  been  achieved.  In  addition  a  nuclear  height  of 
burst,  HOB,  calculation  and  several  high  explosive  detonation 
calculations  were  completed. 

The  nuclear  airblast  program  began  with  a  two- 
dimensional  calculation  of  a  1-kiloton  nuclear  explosion  deto¬ 
nated  at  104  feet  above  the  ground  surface  (Appendix  A)  Fortran 
Computer  Code,  FAST2D,  which  embodied  the  FCT  algorithm  JPBFCT 
(a  derivative  of  ETBFCT  (Boris,  1976)  was  used  to  predict  the 
unsteady  flow  field.  The  calculation  predicted  the  transi¬ 
tion  from  regular  to  Mach  reflection  at  a  ground  range 
approximately  equal  to  the  HOB.  Comparisons  to  shock  tube 
wedge  experiments  showed  similarity  in  the  flow  fields 
(Appendix  B) .  When  compared  to  the  best  available  data  agree¬ 
ment  was  found  to  be  within  201  for  the  high-over-pressure 
(regular  reflection)  region  and  within  104  in  the  double 
peaked  (Mach)  region. 

Further  HOB  investigations  lead  to  simulation  of  high 
explosive  calculations.  HE  Blast  wave  phenomena  include 
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reactive  and  two-phase  flows  associated  with  the  motion  of  the 
chemical  explosion  products.  Additionally,  the  existence  of 
two  shocks,  contact  discontinuity,  and  a  rarefaction  region  make 
the  HE  simulation  more  difficult  numerically  and  physically. 
(Appendix  C)  The  FCT  scheme  in  general  represents  an  accurate 
and  flexible  class  of  methods  for  solving  such  nonsteadv  com¬ 
pressible  flow  problems .  However,  in  models  which  treat  all  the 
physical  effects  required  for  blast  wave  simulation,  truncation 
errors  inherent  in  the  underlying  finite-difference  scheme  are 
exacerbated  by  nonlinear  coupling  between  the  fluid  equations 
and  by  the  greater  complexity  of  the  phenomena  being  simulated. 
The  type  of  errors  are  the  "terraces"  which  develop  in  the 
rarefaction  regions.  Elimination  of  these  errors  became 
necessary  to  accurately  simulate  the  high  explosive  flow 
field.  The  solution  was  to  improve  the  short  wavelength 
phase  and  amplitude  properties  of  the  underlying  algorithm. 

Tests  carried  out  on  scalar  advection  of  simple  density 
profiles  by  a  uniform  flow  field  show  that  terracing  does 
not  require  either  diverging  velocities  or  discontinuities 
in  the  profile,  but  appears  typically  (for  v  >  0)  where  the 
first  and  second  derivatives  of  density  have  the  same  sign. 

In  order  to  improve  the  properties  of  the  basic  difference 
scheme,  we  developed  a  new  algorithm  for  integrating  generalized 
continuity  equations  over  a  timestep  6t.  Consider  the  following 
three-point  transport  scheme: 
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5j  =  C-;  -  T,Co?+1  -  pj.,)  +  <(oJ  +  1  -  2oj  +  Pjc.1); 


^  f  0  0  *v  ^  0  ")0,  o  n 

pj  -  pj  ■  e(pj+i  -  pj-i}  +  x(pi+i  -  2pj  +  pj-i); 


where 


3j  “  pj  vC*j  +  l/2  '  ^ j  - 1/2^ 


*\j  *\* 

3 j  +  1/2  =  pj+l  *  pj- 


The  arrays  {p?}  and  {p*?}  are  the  old  and  new  densities,  Pj 
and  Pj  are  temporary  intermediate  densities,  and  n ,  6,  k, 

A,  and  y  are  velocity-dependent  coefficients.  Here  k  and 
X  are  diffusion  coefficients,  and  y  is  the  antidiffusion 
coefficient.  In  the  actual  algorithm,  4,j+1/2  *s  corrected 

(hence  the  name  FCT)  to  a  value  <^+1/2  chosen  so  no  extrema 

n  ‘  n 

in  Pj  can  be  enhanced  or  new  ones  introduced  in  pj .  Previous 

FCT  algorithms  had  0=0;  the  widely  used  ETBFCT  and  related 

algorithms  (Boris,  1976)  have  in  addition  k  =  0.  If  we 

define  Pj  to  be  sinusoidal  with  wave  number  k  on  mesh  with 

uniform  spacing  5x,  so  that  p?  =  exp  (ijB)  where  $  =  k6x, 

then  the  new  density  array  satisfies 


Pj/Pj  =  A  *  1 - 2 i (n  +  6)sinB  +2  (k  +  X)  (cosg-1) 

-  2y(cos6-l)  [l-2insinB+2ic(cos8-l)]  . 
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From  A  we  can  determine  the  amplification  a  =  A  and  relative 
phase  error  R  =  (l/eS)tan  ^ ( - ImA/ReA) - 1 ,  where  z  =  vct/:x  is 
the  Courant  number.  Expanding  in  powers  of  6  we  find 

a  =  1  +  a76“  +  +  a^B^  +  ...  ; 

R  ■  R0  *  V2  ♦  M4  *  V6  *  •••  ■ 

First-order  accuracy  entails  making  RQ  vanish,  which  requires 
that  n  +  6  =  e/2.  Second-order  accuracy  (a,  =  0)  implies  that 
y  *  k  +  X  -  t~/2.  Analogously,  the  "reduced-phase-error" 
property  R^  =  0  (Boris  and  Book,  1976)  determines  p  =  (l-e")/6, 
thus  leaving  two  free  parameters.  One  of  these  can  be  used  to 
make  R4  vanish  also.  The  resulting  phase  error  R(B)  is  small 
not  only  as  B  -*■  0 ,  but  also  for  larger  values  of  6,  correspond¬ 
ing  to  the  short  wavelengths  responsible  for  terraces.  The 
remaining  parameter  n  can  be  chosen  to  relax  the  Courant  number 
restriction  needed  to  ensure  positivity  from  e  <  1/2  to  t  <  1. 
When  coded,  these  changes  necessitate  a  small  increase  in  the 
operation  count  of  ETBFCT  along  with  a  small  increase  in  over¬ 
head  to  precalculate  the  two  new  arrays  of  velocity-dependent 
transport  coefficients. 

The  technique  described  here  improved  the  high  explosive 
simulation  considerably.  The  decrease  in  phase  error  reduced 
terracing  dramatically  (but  did  not  completely  eliminate  it 
everywhere) .  Smoothing  the  velocity  used  in  the  advection 
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coefficients  and  the  pressure  used  in  the  driving  terms  used 
in  the  advection  coefficients  and  the  pressure  used  in  the 
driving  terms  completely  eliminated  the  problem. 

Next,  a  two-dimensional  (2D)  numerical  calculation  was 
performed  to  simulate  one  of  Carpenter's  (1975)  height-of- 
burst  experiments  which  used  spherical  8-lb.  charges  of 
PBX  9404.  The  previous  fine-zone  ID  calculation  was  used  to 
initialize  the  problem.  It  was  mapped  onto  the  2D  grid  just 
prior  to  the  onset  of  reflection.  The  solution  was  then 
advanced  in  time,  with  pressure  being  calculated  from  a  real- 
air  equation  of  state  and  a  JWL  equation  of  state  for  the 
combustion  products.  The  front  of  the  blast  wave  was  captured 
in  a  finely  gridded  region  which  moved  outward  horizontally. 
Special  care  was  taken  to  ensure  that  the  grid  moved  smoothly. 
The  resulting  solution,  particularly  the  curve  of  peak  over¬ 
pressure  vs.  range,  was  in  excellent  agreement  with  Carpenter's 
experimental  data.  We  believe  that  this  calculation  represents 
the  most  accurate  simulation  of  the  double-Mach-stem  region 
that  occurs  when  a  spherical  blast  wave  reflects  from  a  planar 
surface  which  has  thus  far  been  carried  out. 
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TRANSITION  TO  DOUBLE  MACH  STEM  FOR  NUCLEAR  EXPLOSION  AT 
104  FT  HEIGHT  OF  BURST 


I .  INTRODUCTION 

In  a  nuclear  height-of-burst  (HOB)  detonation  the  spherical  blast 
wave  reflects  from  the  ground,  initially  producing  a  regular  reflection 
region.  When  the  shock  reaches  a  ground  range  approximately  equal  to  the 
HOB  an  abrupt  transition  to  Mach  reflection  occurs.  This  transition  is 
responsible  for  an  airblast  environment  more  severe  than  the  surface  burst 
nuclear  case.  Qualitatively,  it  can  be  thought  of  as  a  partial  flow  stagnation 
in  the  Mach  region  that  leads  to  the  production  of  two  static  pressure  peaks. 

A  1  Kilo ton  (1  KT)  atmospheric  nuclear  explosion  at  a  HOB  of  104  feet  has 
been  simulated  using  the  two  dimensional  FAST2D  code  (Ref  .1).  Figure  1 
illustrates  the  shock  structure.  The  calculation  predicts  the  transition  of  the 
shock  from  regular  reflection  to  double  Mach  reflection.  Because  the 
spherical  waves  are  expanding  and  thus  decreasing  in  Mach  number  as  well  as  angle  of 
incidence  with  the  ground,  they  create  a  dynamic  Mach  stem  formation.  In  comparision 
to  planar  shocks  on  wedges  one  finds  them  to  be  qualitatively  alike.  The  appearance  of 
double  peaks  In  the  pressure  and  density  profiles  (versus  time  and  distance)  is 
interpreted  as  the  point  of  transition.  Other  Interesting  phenomena  such  as  the 
rollup  of  the  contact  surface  generating  a  vortex  ring  and  the  associated 
phenomenon  of  toeing  out  of  the  first  Mach  stem  can  be  observed. 

The  ability  of  the  calculation  to  accurately  predict  the  gasdynamic  effects 
both  temporally  and  spatially  is  due  in  part  to  the  shock  capturing  and  adaptive 
rezone  features  of  the  FAST2D  code.  A  minimal  number  of  very  fine  zones  was  placed 
around  the  shock  front  and  these  zones  then  moved  with  the  first  Mach  stem  to 
prevent  shock  smearing  and  distortion.  This  calculation  is  the  first  sttnpt  to 
model  the  nuclear  HOB  case  through  the  use  of  a  Flux-Corrected  Transport  (FCT) 
algorithm  (Ref.  2). 


Manuscript  submitted  Augwt  24, 1981. 


1 


The  results  of  this  calculation  agree  well  with  the  pressure-distance 
curve  generated  by  the  high  explosive  (HE)  data  of  Carpenter  (Ref.  3j  and  tne 
analysis  of  Kunl  (Ref.  4).  The  peak  overpressure  of  the  first  shock  at  tne  time 
of  transition  is  about  4300  psi.  Our  simulation  was  run  to  11.6  ms  (total 
time  with  t0  -  3.76  ms),  which  corresponds  to  pressure  peaks  of  about  2000  psi. 
In  the  regular  reflection  region  the  peak  values  tend  to  be  about  20*  low  due  tc 
the  clipping  of  the  FCT  algorithm  and  inaccuracies  in  the  initialization  of  the 
flow.  Reducing  the  minimum  zone  size  from  5  cm  to  1  cm  in  a  one-dimensional 
test  calculation  eliminated  this  discrepancy,  however.  In  the  two-peak  regions 
the  agreement  between  the  experimental  data  and  the  values  presented  here  is 
very  good.  The  resolution  of  the  calculation  is  adequate  for  studying  quali¬ 
tatively  the  characteristics  of  the  flow  field.  For  future  work  we  recommend 
that  the  transition  region  be  explored  with  improved  resolution. 
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II.  DESIGN  OF  PROBLEM 


The  problem  of  a  1— K.T  nuclear  detonation  at  104  ft  (31. 7m)  ROB  was  chosen 

since  it  can  be  scaled  conveniently  to  various  HE  tests.  The  use  of  the  1KT 

standard  is  also  expendient;  one  could,  however,  have  used  realistic  initial 

conditions,  such  as  the  Los  Alamos  Scientific  Laboratory  RAD FLO  or  Air  Force 

Weapons  Laboratory  (AFWL)  SPUTTER  calculations.  A  simple  constant  ambient 

-3  3 

atmosphere  was  used  with  a  density  of  1.22  x  10  g/cm  and  pressure  1.01  x 
6  2 

10  dyne/cm  .  To  relate  the  energy  and  mass  densities  to  the  pressure, 
a  real-air  equation  of  state  (EOS)  was  used.  This  "table-lookup"  EOS  is  derived 
from  Gilmore's  data  (Ref.  S.)  and  has  been  vectorized  for  the  TI  Advanced 
Scientific  Computer  at  NRL  (Ref.  6).  Figure  2  illustrates  the  effective  gamma 
versus  specific  energy  per  unit  mass  for  different  values  of  the  density.  The 
internal  energy  density  used  in  the  call  to  the  EOS  is  found  by  subtracting 
kinetic  energy  from  the  total  energy;  this  can  be  negative  due  to  phase  errors 
in  the  fluid  variables.  When  this  occurs,  the  value  of  the  pressure  is  reset  to 
zero. 

The  transition  from  regular  reflection  to  double  Mach  reflection  is  known  to 
occur  at  a  ground  range  approximately  equal  to  the  HOB.  Therefore,  the  size  of  the 
mesh  should  be  roughly  twice  the  HOB  in  both  directions.  The  upper  boundary 

should  be  far  enough  away  from  the  blast  front  to  be  noninterfering.  We  set  the 

3  4 

boundaries  at  3.5  x  10  cm  for  the  radial  direction  and  1.035  x  10  cm  for  the 

axial  direction.  The  fine  grid  in  the  radial  direction  contained  140  out  of 

200  total  zones,  each  5  cm  in  length.  The  largest  zones  initially  filled  the 

right  section  of  the  grid  and  were  80  cm  in  length.  A  smoothing  involving  40 

zones  was  performed  between  the  region  to  guarantee  that  the  zone  sizes  varied 

slowly.  In  che  vertical  direction  the  fine  grid  contained  75  out  of  150  total 

zones,  each  5cm  in  length.  Beyond  that  region  the  zones  increased  geometrically 

by  a  factor  of  1.112. 
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Fig.  2  -  Ganna  -1  vs.  specific  energy  for  Air  EOS 
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Placement  of  the  fine  grid  at  the  origin  (ground  zero  -  the  point  at  which 
first  reflection  occurs)  was  determined  to  be  optimum  for  capturing  peak  pressures 
in  the  airblast  wave  front.  Thus,  as  the  expanding  wave  moved  along  the  ground 
surface,  the  fine  grid  was  always  locked  to  it,  and  each  point  along  the  incident 
blast  front  encountered  the  same  spatial  gridding  as  it  approached  the  ground. 

By  treating  each  point  of  the  incident  front  in  the  same  manner ,  we  insured  that 
the  calculation  was  internally  consistent  and  that  the  computed  transition  point 
was  accurate  to  within  the  limits  of  the  resolution.  Finally,  we  point  out  that, 
as  a  section  of  the  ifet'ient  Mast  wave  propagated  within  the  fine  grid,  the 
wave  steepened.  The  sir  . f  the  fine  grid  was  sufficient  to  insure  that  the  incident 
wave  had  reached  the  maximum  steepeness  prior  to  intersecting  the  ground. 

The  initialization  provides  a  strong  shock  with  Mach  number  M^  ■  12.  This 
speed  and  the  need  for  restart  capability  led  to  the  choice  of  200  timesteps 
as  an  interval  for  the  spatial  display  ("snapshots") .  The  dump  interval  that 
resulted  was  ~At  *  0.3  milliseconds.  These  dumps  were  stored  on  magnetic  tape 
and  postprocessed. 

Additional  diagnostics  were  implemented  in  the  calculation.  Stations  were 
created  to  gain  information  from  fixed  spatial  positions  within  the  calculational 
grid.  These  25  physical  variable  sensors  were  placed  along  the  ground  and  stored 
values  of  the  energy  and  mass  densities  and  velocity  for  every  timestep.  From 
this  information  one  can  construct  static  and  dynamic  pressure  curves. 

III.  COMPUTATIONAL  DETAILS 

The  evolution  of  the  nuclear  HOB  flow  field  was  modeled  numerically  with  the 
FCT  code  FAST2D  (Ref.  7).  FCT  yields  accurate  and  well-resolved  descriptions 
of  shock  wave  propagation  without  the  necessity  of  a  priori  knowledge  of  the 
essential  gasdynamic  discontinuities  in  the  problem.  Additionally,  the  code  has  a 
general  adaptive  regridding  capability  which  permits  fine  zones  to  be 
concentrated  in  the  region  of  greatest  physical  interest  while  the 
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remainder  of  the  system  is  covered  with  coarse  2ones .  Figure  3  depicts  the  grid 
setup  initially  and  at  transition  to  the  double  Mach  stem  structure.  The  rezone 
algorithm  is  programmed  to  track  the  Mach  stem  with  the  fine  grid. 

Tne  transport  algorithm  used  a  low-phase-error  phoenicai  FCT  algorithm  in  a 
model  called  JPBFCT,  an  advanced  version  of  the  ETBFCT  algorithm  described  in  Ref  1. 
The  linear  part  of  this  algorithm  is  fourth-order  accurate  spatially  in  advection 
problems  with  a  giver,  constant  velocity  and  has  a  (nonlinear)  flux-corrected 
antidiffusion  needed  to  model  shocks  correctly.  Finally,  the  transport  subroutine 
is  written  in  sliding-rezone  form,  which  means  that  the  mesh  at  the  beginning  and  the 
end  of  the  timestep  need  not  be  the  same.  Since  the  algorithm  is  one-dimensional, 
timestep  splitting,  is  employed  to  solve  the  2-D  problem. 

The  fluid  transport  routine  JPBFCT  is  fully  vectorized  and  requires  about 
2  us  per  meshpoint-cycle.  This  time  would  have  been  still  less  if  a  vectorized 
fully  two-dimensional  routine  had  been  used,  since  the  1-D  loops  are  too  short  to 
permit  full  advantage  to  be  taken  of  the  vector  capabilities  of  the  NRL  ASC.  The 
table  lookup  in  the  EOS  was  also  fully  vectorized,  so  that  pressure  calculations 
required  about  20Z  of  the  time  needed  for  the  hydrodynamics.  These  two  items 

took  up  nearly  all  of  the  running  time  in  the  blast  calculation  itself.  The  cost 

of  initialization  was  negligible,  but  the  diagnostics  cost  up  to  30X  as  much  as 
the  hydrodynamics,  depending  on  how  many  of  the  various  possible  quantities  were 
actually  plotted.  This  latter  number  would  be  greatly  reduced  if  the  plot  routines 
were  fully  vectorized. 

A  version  of  the  AFWL  1  KT  standard  (Ref  8)  was  used  to  initialize  the  energy, 

density  and  velocity  (flow  field)  at  3.76  milliseconds.  The  corresponding  shock 

4 

radius  was  103.9  ft  (31.69  m)  peak  overpressure  of  1645  psi  (1.134  x  10  K  Pa). 
Because  some  areas  of  the  grid  were  very  coarse,  interpolation  onto  the  grid  was 

performed.  After  the  1  KT  flow  was  laid  down  inside  a  radius  of  104  feet  (31.7  m) 

the  fine-zoned  grid  was  activated  to  follow  the  peak  pressure  as  it  moved  along 
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I 


140  FINE  ZONES  SO  COARSE  ZONES 


Fig.  3  -  Adaptive  gridding.  the  grid  at  initialisation  and  at 
transition  point  (lines  are  drawn  for  every  other  zone,  lines 
in  fine-zone  region  are  indistinguishable). 
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the  ground  surface,  modelled  as  e  perfectly  reflecting  boundary.  This  region 
comprised  HO  zones,  and  a  switch  was  set  to  keep  40  of  these  zones  ahead  of 
the  reflection  point.  Permeable  boundary  conditions  were  used  on  the  top  and 
right  edges  of  the  mesh;  i.e.,  density,  pressure  and  velocity  were  set  equal 
to  ambient  preshock  conditions.  Reflecting  conditions  were  applied  to  the 
left  and  bottom. 

The  timestep  was  recalculated  at  every  cycle  according  to  the  Courant 
condition 

(fac^Ay.) 

dc  ■  0.5  min  TT^fy"!)  ’ 

iJ  (1) 

where  c  is  the  speed  of  sound  and  j V j  is  the  modulus  of  the  flow  speed.  This 
could  have  been  relaxed  somewhat  by  allowing  violation  of  the  local  Courant 
limit  at  points  ("hot  spots")  far  from  the  region  of  chief  physical  interest. 

The  total  elapsed  time  in  the  2-D  calculation,  7.6  ms,  required  5600  cycles. 

Three  types  of  diagnostics  were  employed,  all  in  the  form  of  plots  made 
by  post  nrocessing  a  dump  tape.  The  first  type  of  diagnostic  consisted  of  CRT 
contour  plots  of  density  and  static  pressure,  and  arrows  indicating  the  magnitude 
and  direction  of  the  velocity  field,  obtained  at  the  dump  intervals  (every  200 
cycles).  The  second  type  was  pressure-range  curves  at  z^O,  obtained  by  finding 
the  pressure  peak(s)  along  the  ground  at  each  dump  interval  and  hand-plotting 
them  on  the  same  graph.  The  third  type  consisted  of  pressure  histories  at  a  series 
of  24  stations,  obtained  by  saving  the  energy  and  mass  densities  and  the  veloci¬ 
ties  at  every  cycle. 

IV.  RESULTS  AND  PHENOMENOLOGY 

This  calculation  has  been  done  to  understand  the  violent  effects  of  1  KT 
of  energy  being  released  in  the  atmosphere  at  a  HOB  of  104  ft  (31.7m).  A  strong 
spherical  shock  is  created  in  the  surrounding  air,  and  reflects  from  the  ground. 
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The  outward- traveling  airblast  is  then  composed  of  two  parts:  one  reflected 
upward  approximately  normal  to  the  ground,  and  the  original  spherical  blast. 

The  peak  pressure  is  coincident  with  the  intersection  of  the  two  waves.  This 
intersection  continues  to  move  outward  until  the  angle  of  the  spherical  shock 
with  respect  to  the  ground  reaches  a  critical  value  and  the  transition  to  a 
double  Mach  stem  occurs.  As  shown  by  Ben-Dor  and  Glass  (Ref.  9),  this  angle 
depends  upon  the  incident  strength  of  the  shock  (Fig.  1).  Shocks  with  Mach 
numbers  greater  than  10  are  not  shown.  Initially,  the  Mach  number  for  the  HOB 
simulation  is  well  above  10.  The  Mach  number  at  transition  is  approximately 
11  and  the  angle  is  less  than  50°.  From  Fig.  4  the  corresponding  region  is 
double  Mach  reflection. 

Figure  1  has  been  labeled  with  the  notation  of  Ben-Dor  and  Glass  (Ref.  9). 

It  should  be  noted  that  what  is  generally  regarded  as  the  second  Mach  stem 
is  in  fact  the  second  reflected  wave,  which  is  part  of  the  second  Mach  struc¬ 
ture.  To  be  consistent,  one  must  label  the  second  Mach  stem  as  at  the 

indicated  location.  The  definition  used  is  the  state  of  the  fluid  one 
obtains  by  passing  through  one  shock  wave  (M^)  or  two  shock  waves  (R  and  R^) . 

The  first  reflected  wave  R  becomes  the  incident  wave  for  the  second  Mach 

structure.  Density  contours  are  shown  in  Fig.  5  for  an  planar  shock  on 

wedge  with  a  Mach  number  of  7  and  an  angle  of  50°.  The  complimentary 

figure  illustrates  the  proper  labeling  of  the  multiple  waves.  Comparison 

of  Fig.  5  and  the  HOB  simulation  (Fig.  6)  shows  that  corresponding  waves 

can  be  identified.  Differences  between  the  planar  shock  on  wedge  and  the 

HOB  can  be  explained  in  terms  of  the  unsteady  nature  of  the  HOB  case  (a 

spherically  expanding  wave  that  continuously  decreases  in  Mach  number  and 

angle.)  Although  the  term  irregular  Mach  reflection  has  been  used  to 

describe  the  complex  shock  structure  that  evolves  from  HOB  events,  we  believe  it 

to  be  verv  regular  and  explainable  as  a  double  Mach  reflection  that  evolves  as  a 

function  of  time. 
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The  HOB  numerical  simulation  begins  just  before  the  shock  first  reflects 
from  the  ground.  As  a  summary  of  how  the  flow  field  then  develops,  we  pre¬ 
sent  snap-shots  at  the  important  stages .  (A  more  complete  display  is  presented 
in  Appendix  A) .  Figure  6a  indicates  the  pressure  and  density  contours 
and  velocity  vectors  at  t  «  3,18  ms.  In  Fig.  6b  the  reflected  shock  is  shown 
moving  upward,  the  outward  flow  begins  to  stagnate  at  the  ground  (transition). 

Fig.  6c,  at  t  =  5.99  ms,  shows  an  enlargement  of  the  shock  front,  and  the  development 
of  the  Mach  stem,  slip  surface  and  second  Mach  stem.  The  angle  of  the  shock 
with  respect  to  the  ground  is  increasing  with  time,  so  that  the  effective  wedge 
angle  is  decreasing.  From  the  work  of  Ben-Dor  and  Glass  one  expects  a  transition 
to  double  Mach  stem  to  occur  at  approximately  45°.  The  angle  in  Fig.  6b  is  about 
45°  and  the  shock  front  has  entered  the  transition  phase.  Figure  6d  shows  the 
fully  developed  shock  structure  at  7.79  ms.  Toeing  out  of  the  first  Mach  stem 
can  be  also  seen  in  the  contours  of  Fig.  6d  and  occurs  as  the  fluid  rolls 
forward  where  the  slip  line  would  otherwise  intersect  the  ground.  The  velocity- 
field  in  Fig.  6d  also  shows  this  detail. 

Note  the  reflected  shock  properties  (that  part  of  the  structure  that  contains 
the  second  Mach  stem  M^).  The  reflected  shock  propagates  rapidly  through  the  high 
temperature  fireball,  due  to  the  high  local  sound  speed.  The  shape  of  this  reflected 
wave  is  a  primary  difference  between  the  HOB  case  and  the  planar  wave  on  wedge  case. 

The  other  major  difference,  of  course,  is  the  spherically  expanding  blast  wave  which 

-2 

decreases  in  strength  approximately  as  r 
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PRESSURE 


density 


VELOCITY 


<W|9’"0*  1  *10'*  V*,,.«T«XI05 


Fig.  6  -  Pressure,  density,  and  velocity  fields  for  HOB  calculation 
(a)  in  regular  reflection  state;  (b)  at  transition  to  Mach  reflection; 
(c)  shortly  afterward,  when  second  peak  has  become  larger  than  first; 
and  (d)  fully  developed  (note  toe  at  base  of  first  Mach  stem).  Units 
in  cgs. 
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Ar-  elective  way  to  quantitative-v  evaluate  the  calculation  anc  observe 
ir.  detail  the  transition  to  a  Mach  stem  regime  can  be  seer.  fav  examining  tne 
station  data.  The  station  sensors  were  placed  ir.  the  bottom  rov  c:  tne  cal¬ 
culation  100  to  200  cm  apart.  In  Table  1  the  maximum  pressure  recorded  for 
each  station  along  with  the  location  can  be  found. 

Besides  giving  a  reliable  value  for  the  peak  pressure  to  be  used  for 
constructing  the  pressure-distance  curve,  these  data  allow  one  tc  see  effects 
fixed  in  space  but  varying  in  time.  Figure  7  is  a  superposition  of  the  pressure 
profiles  from  stations  15  to  24  (the  transition  region).  Noteworthy  is  the 
profile  from  station  17,  which  is  the  first  station  to  record  a  second  peak  on 
the  back  side.  At  station  19,  200  cm  further  away  from  ground  aero,  the  second 
peak  is  almost  equal  to  the  first.  The  visible  transition  (as  seen 
in  Fig.  Cb,  occurs  at  a  ground  range  between  3100  cm  and  3300  cm  (stations 
19-21)  revealing  a  dominant  second  peak  and  a  "first  peak"  (i.e.,  first 
seen  by  the  sensor)  that  is  about  half  the  magnitude  of  the  second.  The 
second  peak  does  not  exhibit  a  sharp  almost  discontinuous  rise  and  then  a 
rapid  but  slower  decrease  along  the  back  side.  Instead,  it  has  the  appearance 
of  a  density  compression.  This  behavior  has  dramatic  consequences  for  mi¬ 
litary  planners  because  the  pressure-distance  curve  is  modified  and  the 
dynamic  pressure  is  enhanced. 

The  analogous  profiles  for  dynamic  pressure  are  presented  in  Fig.  8. 

Again  data  from  stations  15  to  24  is  utilized.  The  development  of  the 
second  peak  and  its  correlation  with  the  Mach  stem  formation  can  be  observed. 
There  is,  in  addition,  a  noticeable  increase  in  the  first  peak  values  (station 
15  to  the  maximum  at  station  18).  After  the  structure  becomes  visibly 
resolved  (station  20  and  beyond)  the  second  peak  resembles  a  rounded  profile 
suggesting  the  formation  of  a  stagnation  region  behind  the  first  peak  (Mach 
stem) . 
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Table  1 


Scacion  No. 

Location  (cm) 

Time  (sec^ 

2.0000E  02 

2.25E-04 

2 

4.0000E  02 

2 . 80E-04 

3 

8.0000E  02 

5.28E-04 

4 

1.0000E  03 

7.23E-04 

5 

1.2000E  03 

9.54E-04 

6 

1.4000E  03 

1.23E-03 

7 

1.6000E  03 

1.54E-03 

8 

1.8000E  03 

1 . 91E-03 

9 

2 . 0000 E  03 

2.34E-03 

10 

2.2000E  03 

2.81E-03 

11 

2.3000E  03 

3.07E-03 

12 

2.4000E  03 

3.35E-03 

13 

2.5000E  03 

3.62E-03 

14 

2.6000E  03 

3.88E-03 

15 

2.7000E  03 

4.19E-03 

16 

2.8000E  03 

4.48E-03 

17 

2.9000E  03 

4.79E-03 

18 

3.0000E  03 

5.11E-03 

19 

3.1000E  03 

5.41E-03 

20 

3.2000E  03 

6.06E-03 

21 

3.3000E  03 

6.49E-03 

22 

3.4000E  03 

6.82E-03 

23 

3.5000E  03 

7.28E-03 

24 

3.6000E  03 

7  .73E-03 

Pres  (dynes/cm“) 


8. 11E  08 
7.92E  08 
7.17E  08 
6.73E  08 
6.24E  08 
5.65E  08 
5.21E  08 
4.70E  08 
4.54E  08 
4.14E  08 
4.03E  08 
3.92E  08 
3.82E  08 
3.73E  08 
3.37E  08 
3.33E  08 
3.05E  08 
3.01E  08 
2.33E  08 
1.96E  08 
1.79E  08 
1.87E  08 
1.85E  08 
1.69E  08 


16 


Fig.  8  -  Station  dynamic  preaaura  data  (noa.  15-24) 


Finally  we  consider  the  pressure-distance  relation  for  the  HOE  case. 

In  Fig.  9  we  compare  the  results  of  the  numerical  simulation  with  the  data 
of  Carpenter  and  with  empirical  analysis.  Carpenter's  data  are  based  upon 
careful  HOE  experiments  with  8  lb  PBX9404  spheres.  The  empirical  analvsis 
was  based  on  a  1  KT  nuclear  free  air  curve  and  HOB  construction  factors. 

The  calculated  values  in  the  regular  reflection  regime  are  201  low,  which 
may  be  attributed  to  a  combination  of  FCT  clipping,  the  resolution  of  the  grid, 
and  inaccuracies  in  the  initialization  of  the  flow  field.  During  and  after  Mach 
reflection,  the  peaks  remain  low  until  the  Mach  stem  structure  has  grown  large 
enough  to  be  resolved  on  the  mesh.  By  the  time  it  occupies  a  region  of  15  cells 
high  and  35  cells  wide,  the  peak  pressures  are  in  good  agreement  with  the  HE 
data  and  the  empirical  analysis. 

Other  attempts  to  model  the  transition  region  have  been  made.  Needham  and 
Booen  (Ref.  10)  present  results  of  a  1100  lb  pentolite  sphere  at  15  feet  hOB.  The 
general  phenomena  of  the  flow  field  can  be  seen  from  their  simulation.  When  a 
pressure  distance  curve  is  constructed  from  this  calculation,  one  finds  tnat  in  the 
regular  reflection  region  their  results  are  151  to  301  high  relative  to  theory. 
After  transition  to  double  Mach  reflection  the  first  peaks  are  201  low  while 
the  second  peaks  are  401  low  (Ref.  4). 


V.  SUMMARY 


The  airblast  from  a  1KT  nuclear  event  at  104  ft  HOB  has  been  numerically 
simulated  with  the  FAST2D  computer  code.  The  results  give  insight  to  the 
formation  and  subsequent  evolution  of  the  Mach  stem,  the  triple  point,  and 
the  contact  discontinuity.  The  transition  from  regular  reflection  to  double 
Mach  reflection  is  predicted.  We  suggest  that  the  first  signal  for 
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2D  FCT  CALCULATION  1  KT  AT  104  FT  HOB 

O  FIRST  PEAK 

□  SECOND  PEAK 

THEORETICAL  CURVE:  1  KT  AT  104  FT  HOB 

—A—  1  KT  FREE  AIR  CURVE  ♦  HOB  CONSTRUCTION  FACTORS 
— V—  NUCLEAR  HOB  CURVE  (CARPENTER) 

P B X -  9404 DATA  -CARPENTER  (SHOB  *1.059  FT/LB1/3) 

•  DATA  POINTS  (SHOTS  17,18,19) 


J _ I _ I - L 

10  20  30  40 

GROUND  RANGE  (m) 


9  -  Pressure -range  curves  for  first  and  second  (after  transition 
denoted  by  TP  -  to  double  Mach  reflection)  peaks 


transition  is  the  appearance  of  a  second  peak  behind  the  shock  front  due  to 
stagnation  in  the  flow.  Comparison  with  the  pressure-distance  curves  of 
Carpenter  and  Kuhi  indicates  agreement  within  20%.  Both  first  and  second 
peaks  are  predicted  with  similar  accuracy. 


ACKNOWLEDGEMENT 

We  are  grateful  to  Dr.  Allen  Kuhl  of  R&D  Associates  for  his  valuable 
advice  and  technical  guidance. 

This  work  was  supported  by  the  Defense  Nulcear  Agency  under  Subtask 
Y99QAXSG,  work  unit  00001,  work  unit  title,  "Flux-Corrected  Transport." 


21 


References 


1.  Boris,  J.P.,  Flux-Convected  Transport  Algorithms  for  Solving 
Generalized  Continuity  Equations,  NRL  Report  #3237  (1976) 

2.  Boris,  J.P.,  Book,  D.,  "Flux-Corrected  Transport  I:  SHASTA,  a  Fluid 
Transport  Algorithm  that  Works",  J.  Comp.  Phys.,  11,  38  (1973). 

3.  Carpenter,  H.J.,  Height-of -Burst  Blast  at  High  Overpressure, 

4th  International  Symposium  on  Military  Applications  of  Blast 
Simulation  (1974). 

4.  Kuhl,  A.  (Private  Communication, 1981) . 

5.  Gilmore,  F.R. ,  "Equilibrium  Composition  and  Thermodynamic  Properties 
of  Air  to  24 , 000°K" ,  RM-1543  (Aug  1955). 

6.  Young,  T.R.  (Private  Consnunica t ion,  1981 ) . 

7.  Book,  D.,  Boris,  J.,  Kuhl,  A.,  Oran,  E.,  Picone,  M.,  and  Zalesak,  S. 
Simulation  of  Complex  Shock  Reflections  from  Wedges  in  Inert  and 
Reactive  Gaseous  Mixtures,  Proceedings  of  Seventh  International 
Conference  on  Numerical  Methods  in  Fluid  Dynamics,  Springer-Verlag, 
Stanford,  1980. 

8.  Needham,  et  al . ,AFWL  1  KT  Nuclear  Airblast  Standard,  AFWL  TR-73-55, 

Air  Force  Weapons  Laboratory  (April  1975). 

9.  Ben-Dor.,  Glass  I. I.,  Domains  and  Boundaries  on  Non-Stationarv 
Oblique  Shock-Wave  Reflections.  I.  Diatomic  Gas,  J.  Fluid  Mechanics, 

Vol.  92,  Pt.  3,  pp.  459-496  (1979). 

10.  Booen,  M.W.,  and  Needham,  C.E.,  Two  Dimensional  Hull  Code  Simulation 

of  Complex  and  Double  Mach  Reflections.  Technical  Note  No.  NTE-TN-81-001 , 
Air  Force  Weapons  Laboratory,  New  Mexico,  January  1981. 


APPENDIX  A.  DETAILED  TIME  HISTORY  OF  CALCULATION 


The  following  figures  comprise  e  temporal  history  of  the  numerical 
simulation.  Each  page  contains  pressure  contours,  velocity  vectors, 
density  contours,  and  the  corresponding  grid  for  a  particular  time. 

The  series  begins  at  t  *  0  (tj  *  3.76  ms)  and  continues  to  t^.  «  8.26  ms. 
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A  numerical  technique  has  been  developed  for  capturing  com¬ 
plex,  nonsteady  shock  structures  in  multidimensions.  The 
technique  relies  on  moving  the  computational  mesh  with  the 
shock  wave  so  that  the  features  of  principal  interest  appear 
approximately  stationary.  The  method  has  been  implemented 
using  coordinate-split  Flux-Corrected  Transport  (FCT)  al¬ 
gorithms  which  allow  the  mesh  to  evolve  arbitrarily  with 
respect  to  the  fluid  in  each  coordinate.  The  grid  may  thus 
be  optimized  in  response  to  the  needs  of  a  given  problem. 
Synchronizing  the  grid  and  fluid  motions  permits  significant 
reduction  of  numerical  transients  and  eliminates  numerical 
diffusion.  Shocks  develop  naturally,  with  no  fitting.  The 
method  is  illustrated  by  calculating  complex,  two-dimension¬ 
al  Mach  reflection  phenomena  associated  with  airblasts  and 
shock  diffraction  on  wedges.  The  numerical  results  are  in 
good  agreement  with  available  experimental  data. 

INTRODUCTION 


Numerical  solution  of  transient  multidimensional  gas  dynamics  problems  is 
always  nontrivial.  When,  in  addition,  the  problem  involves  reflecting  super¬ 
sonic  flows,  large  variations  in  length  scales  in  both  space  and  time,  or  phe¬ 
nomena  for  which  neither  analytic  solutions  nor  detailed  experimental  observa¬ 
tions  are  at  hand,  the  state  of  the  computational  art  is  challenged.  Such  a 
problem  arises  in  calculating  the  oblique  reflection  of  shocks  from  solid  sur¬ 
faces  in  planar  geometries  (e.g.  shock  tube  experiments)  or  axisymsnetric  geo¬ 
metries  (e.g.  airblasts).  The  complications  arise  mainly  from  the  presence 
of  Mach  reflections  which  occur  when  a  shock  front  impinges  on  a  reflecting 
surface  at  angles  of  incidence  sufficiently  far  from  normal.  The  formation  of 
a  Mach  stem  and,  consequently,  of  a  slip  surface  intersecting  the  triple  point 
(the  confluence  of  the  incident,  Mach,  and  reflected  waves)  results  from  the 
requirement  that  the  flow  behind  the  reflected  shock  be  parallel  to  the  re¬ 
flecting  surface,  which  cannot  be  achieved  through  regular  reflection. 


Attempts  to  calculate  the  properties  of  the  flow  in  Mach  reflections  date 
back  at  least  to  von  Neumann^  and  the  research  which  grew  out  of  the  wartime 
explosive  studies^-  .  For  the  simplest  problem,  that  of  a  planar  shock 
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reflecting  from  a  plane  surface,  Jones,  Martin,  and  Thornhill  noted  that  it 
is  possible  to  reduce  the  number  of  independent  variables  to  two  by  transform¬ 
ing  to  the  similarity  variables  x/t,  y/t,  a  device  that  was  also  used  by 
Kutler,  et  al®.  Ben-Dor’  developed  a  theory  which  used  shock  polars  to  explain 
some  of  the  features  of  this  problem,  and  solved  the  system  of  algebraic  equa¬ 
tions  obtained  by  combining  the  jump  conditions  across  the  various  disconti¬ 
nuities  (Courant  and  Friedrichs)®  to  describe  the  flow  in  the  neighborhood  of 
the  triple  point.  To  date,  no  satisfactory  treatment  of  the  complete  flow 
field  has  been  published,  although  some  features  (like  the  shape  of  various 
waveforms)  are  quite  easy  to  model. 

In  connection  with  studies  of  both  chemical  and  nuclear  explosions  there 
have  been  many  attempts  to  model  a  spherical  blast  wave  reflecting  from  the 
ground,  the  so-called  height-of-burst  (HOB)  problem.  The  hydrodynamic  pheno¬ 
mena  in  the  two  cases  are  identical,  although  nonideal  effects  (primarily  ex¬ 
plosive  afterburn  in  the  first  instance  and  radiation  preheating  in  the  second) 
are  different.  Previous  attempts  to  model  two-dimensional  complex  shock  re¬ 
flection  have  suffered  from  restriction  to  describing  part  of  the  system,  the 
use  of  a  special  assumption  like  chat  of  self-similarity,  or  less  than  satis¬ 
factory  agreement  with  experimental  data.^ 

The  calculations  discussed  here  represent  a  step  forward  in  overcoming 
these  difficulties.  They  differ  from  previous  numerical  work  in  incorporating 
two  important  computational  developments:  Flux -Corrected  Transport  (FCT)-'-O  and 
an  adaptive  regridding  procedure,  called  "sliding  rezone", H  which  optimizes 
the  mesh  point  distribution  and  hence  the  resolution  of  surfaces  of  disconti¬ 
nuity. 

FCT  is  a  finite-difference  technique  for  solving  the  fluid  equations  in 
problems  where  sharp  discontinuities  arise  (e.g.  shocks,  slip  surfaces  and 
contact  surfaces) .  It  modifies  the  linear  properties  of  a  second-  (or  higher) 
order  algorithm  by  adding  a  diffusion  term  during  convective  transport,  and 
then  subtracting  it  out  "almost  everywhere"  in  the  antidiffusion  phase  of  each 
time  step.  The  residual  diffusion  is  just  large  enough  to  prevent  dispersive 
ripples  from  arising  at  the  discontinuity,  thus  ensuring  that  all  conserved 
quantities  remain  positive.  FCT  captures  shocks  accurately  over  a  wide  range 
of  parameters.  No  information  about  the  number  or  nature  of  the  surfaces  of 
discontinuity  need  be  provided  prior  to  initiating  the  calculation. 

The  FCT  routine  used  in  the  present  calculations,  called  JPBFCT  (an  ad¬ 
vanced  version  of  ETBFCT)^^  consists  of  a  flexible,  general  transport  module 
which  solves  1-D  fluid  equations  in  Cartesian,  cylindrical,  or  spherical  geo¬ 
metry.  It  provides  a  finite  difference  approximation  to  the  conservation  laws 
of  the  general  form: 


q  C  4>dV  »  -f  (u-u  )  •  dA  +  f  xdA  (]_) 

bt  J  oV(t)  J  CA(t)  ~8  *'  dA(t) 

where  represents  the  mass,  momentum,  energy  or  mass  species  in  cell  aV(t), 
u^  and  u  represent  the  fluid  and  grid  velocities,  respectively,  and  t  repre¬ 
sents  tne  pressure/work  terms.  This  formulation  allows  the  grid  to  slide  with 
respect  to  the  fluid  without  introducing  any  additional  numerical  diffusion. 
Thus,  knowing  where  the  features  of  greatest  interest  are  located,  one  can 
concentrate  fine  zones  where  they  will  resolve  these  features  most  effectively 
as  the  system  evolves  (Fig.  1). 

In  the  next  section  we  describe  the  computational  techniques  used  to  solve 
the  wedge  problem  and  present  the  results  of  four  simulations  carried  out  to 
reproduce  experimental  results  of  Ben-Dor  and  Glass. ^  In  Section  III  we  pre¬ 
sent  a  parallel  discussion  for  a  HOB  calculation.  Finally,  in  Section  IV  we 
summarize  our  conclusions. 


Fig.  1.  Adaptive  grids  for  a)  planar  shocks  on  wedge  (double  Mach  shock  fea¬ 
tures  are  indicated);  b)  and  c)  HOB  problem  initially  and  at  transition  point 
(grid  lines  in  fine-zone  region  are  indistinguishable). 

SHOCK-ON-WEDGE  CALCULATIONS 

The  JPBFCT  algorithm  was  used  in  a  2-D  Cartesian  version  of  the  FAST2D  code 
to  model  the  reflections  of  planar  shocks  from  wedges  of  20°  to  60°  and  vary¬ 
ing  shock  strengths.  Four  general  classes  which  include  regular,  single,  com¬ 
plex  and  double  Mach  reflection  were  calculated  (referred  to  as  cases  a,b,c,d 
respectively).  The  bottom  of  the  mesh,  treated  as  a  reflecting  boundarv, 
modeled  the  surface  of  the  wedge.  Quantities  on  the  right  hand  boundary  and 
on  the  top  were  set  equal  to  the  ambient  values.  The  remaining  boundaries  were 
treated  as  permeable.  In  the  single,  complex,  and  double  Mach  reflection 
cases,  the  mesh  was  anchored  on  the  left,  essentially  at  the  wedge  tip  where 
the  incident  shock  first  strikes,  while  the  zones  were  stretched  by  a  scaling 
factor  proportional  to  t  as  soon  as  the  reflection  region  filled  a  substantial 
portion  of  the  grid.  In  case  (d) ,  the  double  Mach  reflection  case,  the  open¬ 
ing  angle  is  so  small  that  the  incident  shock  has  to  traverse  many  zones 
before  the  mach  stem  has  grown  large  enough  to  be  well  resolved.  For  this  rea¬ 
son,  the  problem  was  solved  on  a  uniform  mesh  in  the  frame  of  reference  fixed 
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to  the  reflection  point,  with  stretching  being  initiated  after  the  first  Mach 
stem  reached  'v-  20  cells  in  length.  The  timestep  was  recalculated  at  every 
cycle  with  a  Courant  number  of  0.5. 

Figure  2  shows  the  pressure  and  density  contours  and  the  velocity  field  for 
cases  a,b,c,d.  The  pertinent  shock  phenomena  can  be  easily  identified:  inci¬ 
dent  shock,  contact  surface,  first  and  second  Mach  stems.  As  shown  in  Fig.  1, 
the  zoning  is  particularly  sparse  except  for  the  region  of  interest.  Adequate 
resolution  of  the  key  surfaces  (contact  and  second  Mach  stem)  is  obtained  with 
5  zones  in  each  direction.  The  accuracy  can  be  evaluated  by  comparing  the  ex¬ 
perimental  density  distributions  along  the  wall  (Fig.  3). 


PRESSURE 


OENSITY 


VELOCITY 


Fig.  2  -  Pressure  and  density  contours  and  flow  velocity  vectors  (in  frame  of 
reflection  point)  for  planar  waves  with  Mach  number  M  reflecting  from  wedges 
with  angle  9  for  (a)  M-2.03,  9-60°;  (b)  M-2.82,  9-20°;  (c)  M-5.29;  9-30°; 

(d)  M-7.03;  9-50° 
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Fig.  3.  Comparison  of  density  (in  units  of  ambient  density  c  )  for  cases  (a), 
(b),  (c),  (d)  of  Fig.  2  vs.  distance  from  comer.  Points  are  measured  values 
reported  in  Ref.  13. 

HEIGHT  OF  BURST  CALCULATIONS 

Next,  we  performed  a  numerical  simulation  of  a  1KT  nuclear  detonation  at 
31.7  m  HOB,  a  case  which  could  be  readily  compared  with  high  explosive  data.  A 
constant  ambient  atmosphere  was  used  with  a  density  of  1.22  x  10“3  g/cm^  and 
pressure  1.01  x  10^  dynes/cm^.  To  relate  the  energy  and  density  to  the  pressure, 
a  real-air  equation  of  state  (EOS)  was  used.  This  table-lookup  EOS  was  derived 
from  theoretical  calculations  by  Gilmore^* ^  for  equilibrium  properties  of  air 
and  has  been  vectorized  for  the  Advanced  Scientific  Computer^  .  The  internal 
energy  density  used  in  the  call  to  the  EOS  is  found  by  subtracting  kinetic 
from  total  energy;  this  can  be  negative  due  to  truncation  (phase)  errors.  When 
this  occurred,  the  value  of  the  pressure  was  reset  to  zero. 

The  transition  from  regular  reflection  to  double  Mach  reflection  occurs 
at  a  ground  range  approximately  equal  to  the  HOB.  The  size  of  the  mesh  should 
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therefore  be  roughly  twice  the  HOB  in  both  directions.  The  upper  boundary 
should  be  tar  enough  away  from  the  blast  front  to  be  non-interfering.  We  chose 
boundaries  of  55  m  for  the  radial  direction  and  103.5  m  for  the  axial  direction. 
The  fine  grid  in  the  radial  direction  contained  140  out  of  200  total  zones, 
each  5  cm  in  length.  The  rightmost  zones  were  80  cm  in  length,  and  a  smoothing 
involving  40  zones  was  performed  between  the  regions  to  guarantee  that  the  zone 
sizes  varied  slowly.  In  the  axial  direction  the  fine  grid  contained  75  out  of 
150  total  zones,  each  5  cm  in  length.  Beyond  that  region  the  zones  were  geo¬ 
metrically  increased  by  a  factor  of  1.112. 

Placement  of  the  fine  grid  at  the  origin  of  the  mesh  (ground  zero,  the 
point  at  which  reflection  first  occurs)  was  determined  to  be  optimum  for  cap¬ 
turing  peak  pressure  in  the  airblast  wavefront.  Thus,  as  the  expanding  wave 
moves  along  the  ground  surface,  the  fine  grid  is  always  locked  to  it  and  each 
point  along  the  blast  front  encounters  the  same  spatial  gridding  as  it  approaches 
the  ground.  8v  treating  each  point  of  the  incident  front  in  the  same  manner, 
we  insure  that  the  calculation  is  internally  consistent  and  that  the  computed 
transition  point  is  accurate  to  within  the  limits  of  the  resolution. 

The  initialization  provides  a  strong  shock  with  approximate  Mach  number 
M=12.  This  speed  and  the  need  for  restart  capability  led  to  the  choice  of  200 
timesteps  as  an  interval  for  the  spatial  display  (snapshots).  The  dump  inter¬ 
val  that  resulted  was  At  0.3  milliseconds  (ms).  These  dumps  were  stored  on 
magnetic  tape  and  post-processed. 

A  fit  to  the  1-D  nuclear  blast  flow  field  (Ref.  17)  was  used  to  initialize 
the  energy  and  mass  density  and  velocity  field  at  3.76  ms.  The  corresponding 
peak  overpressure  was  113  bars.  After  the  1  KT  flow  field  was  laid  down  inside 
a  radius  of  31.6  m,  the  fine-zone  grid  was  activated  to  follow  the  peak  pressure 
as  it  moved  along  the  ground  surface,  modelled  as  a  perfectly  reflecting  bound¬ 
ary.  This  region  comprised  140  zones,  and  a  switch  was  set  to  keep  40  of  these 
zones  ahead  of  the  reflection  point.  Permeable  boundary  conditions  are  used 
on  the  top  and  right  edges  of  the  mesh,  i.e.,  density,  pressure  and  velocity 
are  set  equal  to  ambient  preshock  conditions.  Reflecting  conditions  were 
applied  to  the  left  and  bottom.  The  total  elapsed  physical  time  in  the  2-D 
calculation,  7.6  ms,  required  5600  cycles.  Times  are  referred  to  t=0  at  the 
start  of  the  calculation. 


The  numerical  simulation  begins  just  before  the  shock  first  reflects  from 
the  ground.  Fig.  4a  indicates  the  pressure  and  density  contours  and  velocity 
vectors  at  time  3.18  ms.  In  Fig.  4b  the  reflected  shock  is  shown  moving  upward, 
the  outward  flow  begins  to  stagnate  at  the  ground  (transition).  Fig.  4c,  t*5.99 
ms,  shows  an  enlargement  of  the  shockfront,  and  the  development  of  the  Mach 
stem,  slip  surface  and  second  Mach  stem.  The  angle  of  the  shock  front  with 
respect  to  the  ground  is  increasing  with  time,  so  that  the  effective  wedge 
angle  is  decreasing.  From  Ben-Dor  and  Glass-*-®  one  expects  a  transition  to 
double  Mach  stem  to  occur  at  approximately  45°.  The  angle  in  Fig.  4b  is  about 
45°  and  the  shock  front  has  entered  the  transition  phase.  Figure  4d  shows  the 
fully  developed  shock  structure  at  7.79  ms.  Clearly  visible  is  the  second  Mach 
stem  and  a  vortex  region  behind  the  first  Mach  stem.  Toeing  out  of  the  first 
Mach  stem  can  be  also  seen  in  the  contours  of  Fig.  4d  and  occurs  as  the  fluid 
rolls  forward  where  the  slip  line  would  otherwise  intersect  the  ground.  The 
velocity  field  in  Fig.  4d  also  shows  this  detail. 


One  should  also  note  the  reflected  shock  properties.  The  reflected  shock 
propagates  rapidly  through  the  high  temperature  fireball,  due  to  the  high  local 
soundspeed.  The  shape  of  this  reflected  wave  is  a  primary  difference  between 
the  HOB  case  and  the  wedge  case*-^.  The  other  major  difference,  of  course,  is 
the  spherically  expanding  blast  wave  which  decreases  in  strength  approximately 
proportional  to  r~“. 
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Fig.  4.  Pressure,  density,  and  velocity  fields  for  HOB  calculation  (a)  in 
regular  reflection  stage;  (b)  at  transition  to  Mach  reflection;  (c)  shortly 
afterward,  when  second  peak  has  become  larger  than  first;  and  (d)  fully  de¬ 
veloped  (note  toe  at  base  of  first  Mach  stem) . 

Finally  we  consider  the  pressure/distance  relation  for  the  HOB  case.  In 
Fig.  5  we  compare  the  results  of  the  numerical  simulation  with  the  data  of 
Carpenter  and  with  empirical  analysis.  Carpenter's  data  are  based  upon  care¬ 
ful  HOB  experiments  with  8  lb  PBX9404  spheres.  The  empirical  analysis  was 
based  on  a  1  KT  nuclear  free  air  curve  and  HOB  construction  factors.  The  cal¬ 
culated  values  in  the  regular  reflection  regime  are  20%  low  and  may  be  attri¬ 
buted  to  a  combination  of  FCT  clipping,  the  resolution  of  the  grid,  and  in¬ 
accuracies  in  the  initialization  of  the  flow  field.  During  and  after  Mach 
reflection,  the  peaks  remain  low  until  the  Mach  stem  structure  has  grown  large 
enough  to  be  resolved  on  the  mesh.  By  the  time  it  occupies  a  region  of  15  cells 
high  and  35  cells  wide,  the  peak  pressures  are  in  good  agreement  with  the  HE 
data  and  the  empirical  analysis. 
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GROUND  RANGE  (m) 

Fig.  5.  Pressure-range  curves  for  first  and  second  (after  transition  -  denoted 
by  TP  -  to  double  Mach  reflection)  peaks. 

SUMMARY  AND  CONCLUSION 

The  complex  2-D  Mach  reflection  phenomena  associated  with  shock  diffraction 
on  wedges  and  height-of-burst  explosions  have  been  modeled  with  the  FAST2D  com¬ 
puter  code.  Four  wedge  cases  —  regular,  single,  complex  and  double  Mach  reflec¬ 
tion — have  been  calculated  and  the  results  compared  to  experiments.  A  nuclear 
detonation  (1  KT  at  31.7m  HOB)  was  also  simulated.  The  results  give  insight  in¬ 
to  the  formation  and  subsequent  evolution  of  the  Mach  stem,  the  triple  point 
and  the  contact  discontinuity.  The  transition  from  regular  reflection  to  double 
Mach  reflection  is  predicted.  Excellent  agreement  with  Ben-Dor's  data  is  ootained. 
We  suggest  that  the  first  signal  for  transition  is  the  appearance  of  a  second 
peak  behind  the  shock  front  due  to  stagnation  in  the  flow.  Calculated  first  and 
second  pressure  peaks  versus  distance  in  the  HOB  case  agree  both  with  the  HE 
data  and  analysis  to  within  20%. 

The  use  of  the  adaptive  regridding  procedure,  called  "sliding  rezone", along 
with  the  FCT  algorithm  allows  one  to  accurately  predict  the  nonsteady  shock 
structures  in  two  dimensions  for  diffractions  on  wedges  and  HOB  cases.  Compari¬ 
son  with  data  for  both  wedges  and  HOB  yields  the  best  results  obtained  to  date. 
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20.  Abstract  (Continued) 

contains  multiple  peaks  with  enhanced  loads  and  impulses.  There  is  an  ongoing  interest 
in  simulating  these  HOB  environments  for  military  applications.  High  explosives  (HE) 
charges  can  be  used  to  simulate  the  nuclear  surface  burst  case  below  about  100  psi  for 
reasonable  yields  (100T  or  more),  but  it  appears  that  it  is  impractical  to  elevate  large 
HE  charges  above  grade  to  simulate  the  HOB  case.  In  this  paper  we  propose  a  new 
method  for  naturally  simulating  such  HOB  environments  on  a  large  scale.  A  hemi¬ 
spherical  HE  charge  could  be  detonated  near  a  natural  slope  which  had  been  graded  to 
form  a  large  ramp.  When  the  spherical  blast  wave  reflects  from  this  ramp  a  shock 
structure  and  environment  is  created  which  is  similar  to  the  HOB  case.  Validity  of  this 
concept  is  demonstrated  by  numerical  simulations  with  a  nonsteady  2-D  FCT  hydrocode. 
These  calculations  indicate  that  a  30s  ramp  located  200  ft  from  a  500T  hemispherical 
HE  charge  will  create  400  to  600  psi  double  peak  static  pressure  waveforms  at  distances 
of  40  to  60  ft  up  the  ramp;  time  between  peaks  is  1  ms.  These  waveforms  correspond 
to  a  nuclear  detonation  at  100  to  120  ft/KT1'3  and  a  ground  range  of  190  to  210  ft/ 
KT1'3. 
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FCT  SIMULATION  OF  HOB  AIRBLAST  PHENOMENA 


X .  INTRODUCTION 

It  is  now  recognized  that  height-of-burst  (HOE)  detona¬ 
tions  can  create  more  severe  airblast  environments  than 
surface  burst  (SB)  detonations,  especially  at  high  over¬ 
pressures.  In  the  HOB  case,  the  spherical  blast  wave  re¬ 
flects  from  the  ground,  initially  as  a  regular  reflection. 
Then  at  a  ground  range  approximately  equal  to  the  height- 
of-burst,  the  shock  reflection  makes  a  transition  to  a 
double  Mach  shock  structure.  This  double  shock  structure 
creates  secondary  peaks  in  the  static  pressure  at  and  near 
the  ground  and  thus  enhances  the  early-time  HOB  airblast 
impulses  compared  to  the  SB  case.  As  shown  experimentally 
by  H.  J.  Carpenter  at  MABS-IV  (Ref.  1),  these  secondary 
peaks  of  the  HOB  case  can  be  much  greater  than  the  first 
peaks. 

When  a  double  Mach  shock  structure  reflects  from  an 
above-ground  structure,  it  can  produce  enhanced  diffraction 
loads.  HOB  diffraction  loads  are  compared  with  SB  loads  In 
Fig.  1  which  was  constructed  by  scaling  data  from  the 
1000-lb  Pentolite  sphere  experiments  on  the  recent  MIGHTY 
MACH  test  series  (Ref.  2).  As  is  evident  from  this  figure, 
the  early-time  HOB  loading  impulses  are  about  twice  the 
SB  values.  Similar  effects  are  shown  in  Fig.  1  for  the 
static  pressure  histories  and  impulses  which  apply  to  loads 
on  flush  mounted  structures. 

For  military  applications,  there  is  a  need  to  simulate 
these  HOB  blast  environments  on  8  large  scale  in  order  to 
test  the  response  and  survivability  of  large-scale  or  full- 
scale  military  systems.  Explosive  yields  from  kilotons  to 
megatons  are  required.  Suspension  of  such  large  high  explo¬ 
sive  (HE)  charges  is  impractical  and  could  lead  to  poor 
quality  blast  fields  due  to  interference  effectB  from  the 
charge  support  structure. 

In  this  paper  we  propose  a  novel  approach  for  simulating 
HOB  blast  environments  on  a  large  scale.  The  concept  is 
shown  in  Fig.  2.  A  hemispherical  surface  burst  HE  charge 
would  be  used  to  create  a  free-field  blast  wave.  The  charge 
would  be  situated  near  an  up-slope  which  had  been  graded  to 
form  a  large  ramp.  When  the  spherical  blast  reflects  from 
the  ramp,  a  double  Mach  shock  structure  can  be  created 
(within  certain  constraints  on  wedge  angle,  6y,  and  inci¬ 
dent  shock  Mach  number) .  This  concept  relies  on  the  simi¬ 
larity  between  the  HOB'produced  environments  on  horizontal 
surfaces  and  the  environments  produced  by  shock  reflections 
on  wedges  or  ramps.  In  Fig.  3  we  compare  some  recent 
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calculations  with  the  FAST2D  code:*  a  nuclear  detonation  at 
HOB  »  104  ft/KTjj  ' ^  versus  a  Mach  seven  square  wave  shock 
reflection  from  a  wedge.  The  pressure  contours  show  chat 
for  similar  shock  strengths  and  angles,  the  shock  structures 
in  the  wedge  and  HOB  cases  are  qualitatively  similar;  den¬ 
sity  contours  are  also  qualitatively  similar  with  a  slip 
line  emanating  from  the  primary  triple  point.  There  are, 
however,  quantitative  differences:  the  Mach  stem  structure 
in  the  nuclear  HOB  case  is  more  complex,  with  a  bulge  at 
the  foot  of  the  Mach  stem;  also,  in  the  nuclear  case,  the 
reflected  shock  races  rapidly  through  the  high  temperature 
(lO^  to  10J  K)  fireball,  while  in  the  wedge  case,  the 
reflected  wave  propagates  slowly  into  the  lower  (^lO^  K) 
temperature  constant  field  behind  the  incident  square  wave 
shock.  We  believe,  however,  that  these  differences  are  of 
secondary  importance. 

A  remaining  question  is:  how  well  does  the  blast  wave 
from  a  hemispherical  HE  charge  simulate  the  nuclear  free- 
field  environment?  In  Figure  4  we  compare  the  static  and 
dynamic  pressure  waveforms  for  the  HE  and  nuclear  cases 
from  Brode's  one-dimensional  (l-D)  free  air  burst  calcula¬ 
tions  (Refs.  3,4)  at  shock  overpressures  of  approximately 
100,  200  and  400  psig.  In  the  HE  case  a  contact  surface 
(CS)  separates  the  air  from  the  detonation  products.  This 
contact  surface  causes  a  sharp  jump  in  dynamic  pressure  due 
to  the  high  densities  of  the  products.  Also  evident  in  the 
HE  case  is  a  secondary  shock,  S2,  which  faces  inward  but  is 
being  swept  outward  by  the  .apid  expansion  of  the  charge. 

The  HE-driven  blast  wave  gives  a  rather  poor  simulation  of 
the  complete  nuclear  waveform  at  high  overpressures,  due 
principally  to  the  HE  contact  surface  and  secondary  shock. 
However,  the  HE  blast  wave  outside  the  contact  surface  is  a 
reasonably  good  simulation  of  the  nuclear  case.  We  propose 
to  use  precisely  this  part  of  the  HE  blast  wave  and  reflect 
it  from  the  ramp  to  simulate  the  early-time  nuclear  HOB 
cases. 

The  remainder  of  the  paper  is  organized  as  follows: 
Section  II  gives  a  conceptual  design  of  the  ramp  HOB  simu¬ 
lator;  Section  III  describes  the  2-D  finite  difference 
scheme  which  we  used  to  investigate  numerically  the  flow 
fields  on  and  near  the  ramp;  Section  IV  presents  the  results 
of  these  calculations,  while  conclusions  and  recommenda¬ 
tions  are  offered  in  Sections  V  and  VI. 


* 

This  code  uses  the  Flux  Corrected  Transport  (FCT)  algor¬ 
ithm  .described  in  Section  III,  to  maintain  sharp  disconti¬ 
nuities. 
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II.  CONCEPTUAL  DESIGN  OF  THE  HOB  SIMULATOR 


The  design  objective  for  this  simulator  is  to  produce 
the  high  overpressure  (say  100  to  1000  psi)  double-peak, 
flow  fields  which  simulate  nuclear  HOB  detonations  in  the 
Mach  reflection  regime  with  high  fidelity.  The  simulator 
should  be  reasonably  inexpensive  and  readily  constructed. 
The  design  concept  should  be  extendable  to  large  yields. 

The  primary  design  parameters  for  the  simulator  are  the 
location  of  the  front  edge  of  the  ramp,  GR^,  and  the  ramp 
angle,  9y .  The  conceptual  design  process  begins  with  an  HE 
free-air  pressure-range  curve  for  1  lb  of  Pentolite.  A 
ramp  was  assumed  to  be  located  at  a  GRr  corresponding  to 
free  field  shock  overpressures  of  500  psi  or  150  psi. 
Assuming  various  ramp  angles,  we  used  reflection  factors 
(Ref.  1)  to  determine  the  peak  static  pressure  versus  ramp 
ground  range,  RGR.  Parametric  results  are  presented  in 
Fig.  5.  Inserts  give  the  results  scaled  to  a  500T 
surface  burst  which  are  equivalent  to  about  a  one-kiloton 
nuclear  surface  burst  case.  Examination  of  the  results  in 
Fig.  5  indicates  the  following  trends: 

o  A  requirement  for  a  high  pressure  (400  psi  to  600 
psi)  simulator  forces  one  to  either  move  the  ramp 
closer  to  the  charge,  or  increase  the  ramp  angle, 
or  both. 

o  One  vould  prefer  to  move  the  ramp  away  from  the 
charge  so  that  the  HE  free  field  is  close  to  the 
nuclear  case;  however,  this  leads  to  large  (and 
presumably  impractical)  ramp  angles. 

o  Decreasing  the  ramp  angle  tends  to  make  the  Mach 
stem  rise  more  rapidly  thus  increasing  the  separa¬ 
tion  between  the  first  and  second  peaks;  we  specu¬ 
late  that  this  could  lead  to  a  yield  amplification 
on  the  front-end  of  the  waveform. 

o  Transition  to  Mach  reflection  occurs  at  the  leading 
edge  of  the  ramp  for  9y  -  30  and  40  ;  the  transi¬ 
tion  point  (TP)  for  the  9y  *  60  occurs  at  about 
one-half  the  distance  up  the  ramp. 

The  30°  ramp  at  the  500  psi  station  appears  to  be  an 
interesting  case — it  is  feasible  to  construct  and  the  600- 
psi  shock  overpressure  will  occur  at  about  50  ft  up  the 
ramp, thus  allowing  plenty  of  time  for  the  Mach  stem  height 
to  gTOw.  Peak  pressures  will  range  from  1500  psi  at  the 
beginning  of  the  ramp  to  about  300  psi  at  the  far  end. 
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A  numerical  simulation  of  the  shock  diffraction  for  the 
ramp  HOB  simulator  (a  30  ramp  starting  at  200  feet  from  a 
SOOT  hemispherical  HE  charge)  was  performed  with  a  nons~eady 
two-dimensional  (2-D)  hydrocode,  FAST2D.  The  objectives  of 
the  calculation  were  to  validate  the  ramp  HOB  simulator 
design  and  to  evaluate,  in  detail,  the  flow  field  in  the 
vicinity  of  the  ramp.  The  FAST2D  code  solves  the  balance 
laws  of  gasdynamics  on  a  sliding  grid  in  the  general  form: 


3t 


J  <fdV 
5V(t) 


-e  ^  (u-u 
•^(t)  “8 


dA  +  6  TdA 
y5A(t) 


(1) 


where  $  represents  the  mass,  momentum,  energy  or  species  mass 
density  (for  multi-material  calculations)  in  cell  t5V(t), 
u  and  ug  represent  the  fluid  and  grid  velocities,  respect¬ 
ively,  and  x  represents  the  pressure/work  terms.  The  finite- 
difference  approximation  to  Eq.  (1)  uses  a  vectorized  Flux- 
Corrected  Transport  (FCT)  algorithm,  ETBFCT  (Ref.  5),  which 
gives  an  accurate  and  well-resolved  description  of  shock 
wave  propagation  without  the  necessity  of  an  £  priori  know¬ 
ledge  of  the  number,  location  or  character  of  the  gas- 
dynamic  discontinuities  in  the  problem.  The  linear  portion 
of  this  algorithm  is  fourth-order-accurate  spatially  for 
constant-velocity  advection  problems,  and  has  a  nonlinear 
flux-corrected  antidiffusion  stage  which  automatically  pro¬ 
vides  the  local  dissipation  needed  to  accurately  model  dis¬ 
continuities.  The  formulation  of  the  algorithm  allows  the 
grid  to  slide  with  respect  to  the  fluid  without  introducing 
additional  numerical  diffusion.  This  general  adaptive 
regridding  technique  permits  fine  zones  to  be  concentrated 
in  the  region  of  greatest  physical  interest,  thus  reducing 
computational  costs  with  no  serious  loss  in  resolution. 


Since  the  ETBFCT  algorithm  is  one-dimensional, 
time -splitting  must  be  employed  to  solve  two-dimensional 
problems.  Time- split ting  makes  the  boundary  condition  on 
the  ramp  particularly  easy  to  implement.  The  ramp  is 
represented  as  a  series  of  "stairsteps"  (of  varying  height 
and  depth)  along  the  interface  between  the  extremal  interior 
zones  and  a  corresponding  set  of  guard  cells.  A  guard  cell 
is  defined  as  the  right-most  cell  in  the  r-direction 
during  the  r-sweep,  and  the  bottom-most  cell  in  che  z- 
direction  during  the  z-sweep.  The  stairstep  boundary  con¬ 
ditions  are  ref lective, which  requires  pressure,  density 
and  energy  to  be  continuous  and  the  corresponding  velocity 
normal  to  the  stairstep  to  vanish. 
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The  numerical  simulation  began  with  a  1-D  FCT  calcula¬ 
tion  of  the  blast  wave  driven  by  a  one  pound  spherical 
charge  of  PBX-9404  in  air.  The  initial  conditions,  which 
are  shown  in  Fig.  6,  were  taken  to  be  the  self-similar 
flow  field  corresponding  to  a  spherical  Chapman-Jouguet 
detonation  wave  (Ref.  6),  at  the  time  the  detonation  wave 
reaches  the  charge  radius,  r0  «  3.89  cm/lb^'  .  A  Jones- 
Wilkins-Lee  (JWL)  equation  of  state  (EOS)  was  used  for 
the  detonation  products  and  a  real  air  equation  of  state 
was  used  outside  the  HE/air  interface.  These  EOS  specify 
the  pressure  as  a  function  of  density  and  internal  energy. 
The  HE/air  interface  was  followed  by  solving  a  conserva¬ 
tion  law  for  the  mass  fraction  pf  (where  f*l  in  the  pure  HE 
and  f*»0  in  the  pure  air).  The  equations  of  state  were 
blended  in  the  mixed  cells  (0<f<l)  according  to  Dalton's 
law.  A  fixed  grid  of  500  cells  was  used  with  a  mesh  spacing 
Ar  *  0.1025  cm/lb3'3,  so  that  the  initial  flow  field  in  the 
charge  occupied  about  38  computational  cells.  The  flow 
field  results  at  the  end  of  the  1-D  calculation  (cycle  1281, 
t  *  152  ps/lb3'3)  are  shown  in  Fig.  6.  The  shock  over¬ 
pressure  is  445  psig.  The  density  distribution  shows  a 
jump  at  the  HE/air  interface;  inside  the  interface  is  a 
secondary  inward-facing  shock  which  is  being  swept  outward 
by  the  supersonic  flow. 

These  results  were  scaled  up  to  the  500  ton  HE  surface 
burst  case  by  multiplying  all  times  and  ranges  by  the  scale 
factor,  SF  »  (2x10^) 1/3  «  125.992.  The  shock  radius  at 
this  time  of  19.15  ms/SOOT1/3  was  found  to  be  198  £t/500T1/3 
with  an  overpressure  of  445  psig  (note  that  this  point 
checks  with  the  HE  free  air  curve  in  Fig.  5).  These 
results  were  then  inserted  as  initial  conditions  in  the 
cylindrical  r-z  FAST2D  code,  with  one  approximation.  Since 
the  y’s  ahead  and  behind  the  HE/air  interface  were  quite 
close  (yhe  *  1.25  versus  Ya^r  *  1.30),  the  HE  products  were 
modeled  with  the  real  air  equation  of  state,  and  the  2-D 
interface  was  not  followed  specifically  with  a  mass  species 
conservation  law.  The  2-D  mesh  consisted  of  150  x  150 
cells  with  a  moving  fine  mesh  region  (55  x  55  cells,  Ar  ■ 

5  cm  and  Az  *■  2.8868  cm  with  Az/Ar  »  tan  30°)  which  followed 
the  Mach  stem.  The  calculation  was  run  5601  cycles.  Diag¬ 
nostics  for  the  2-D  calculation  consisted  of  46  environment 
time  histories  (at  40  stations  on  the  ramp  and  6  stations 
perpendicular  to  the  ramp  at  a  RGR  -  60.5  ft)  and  contour 
plots  of  the  flow  field  every  200  cycles.  Times  are 
denoted  by  the  label  At-t-c0j  which  references  everything 
to  the  incident  shock  arrival  time  at  the  foot  of  the  ramp 
t0*19.2  ms. 
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IV.  CALCULATIONAL  RESULTS 


An  overall  picture  of  the  spherical  shock,  reflection 
from  the  ramp  is  displayed  in  Fig.  7  which  gives  the  cal¬ 
culated  pressure  and  density  contours  at  various  times 
(!t=3 ,  5.61,  9.27  and  13.4  ms/SOOT1-'3)  ;  Fig.  8  gives  a 
magnified  view  of  the  flow  field  at  2t=9.61  ms/500T3''3. 

The  shape  of  the  shock  structure  for  the  simulator  ("i.e., 
the  geometry  of  the  incident  wave,  the  Mach  stem,  and  the 
kinked  reflected  wave)  more  closely  resemble  the  shock 
structure  for  square  wave  reflections  from  a  ramp  (Ref.  7) 
than  the  nuclear  HOB  case  (see  Fig.  3).  The  density 
contours  indicate  that  a  contact  surface  (a  slip  line) 
emanates  from  the  triple  point  (the  confluence  of  the  inci¬ 
dent,  Mach  and  reflected  waves)  and  approaches  the  ramp  at 
an  angle  of  about  60  .  Pressure  contours  indicate  that  a 
high  pressure  region  is  located  in  the  vicinity  of  where 
the  projection  of  the  contact  surface  would  strike  the  ramp. 

Figure  9  gives  an  experimental  shadowgraph  of  the  shock 
wave  structure  formed  by  an  3  lb  TNT  driven  blast  wave 
(HOB  »  1.04  ft/lb3-/3)  diffracting  on  a  31°  ramp.  The  inci¬ 
dent  shock  pressure  was  about  120  psi  at  the  foot  of  the 
ramp  and  about  75  psi  at  the  time  of  the  photograph  (com¬ 
pliments  of  W.  Dudziak,  Ref.  8).  The  shock  structure  is 
qualitatively  similar  to  that  in  Figs.  7  and  3.  Fig.  9 
shows  that  the  reflected  wave  pushes  the  TNT  products  away 
from  the  ramp,  thus  maintaining  a  clean  air  flow  (unpolluted 
by  HE  products)  in  the  Mach  stem  region — a  truly  beneficial 
result!  Note  that  this  happens  even  in  the  low  HOB  case 
where  the  TNT  products  squish  along  the  ground  and  push  the 
TNT/air  interface  closer  to  the  shock. 

The  calculated  shock  properties  for  the  ramp  HOB  simu¬ 
lator  are  shown  in  Fig.  10  as  a  function  of  ramp  ground 
range,  RGR.  The  primary  Mach  stem  pressure,  pj_,  ranged 
from  about  600  psi  to  400  psi.  The  second  peak  pressure, 

P2,  decayed  from  1300  psi  at  the  foot  of  the  ramp  to  400 
psi  at  the  60  foot  station.  The  peak  pressures  were  deter¬ 
mined  from  two  methods:  for  RGR  <  30  ft  peaks  were  evaluated 
from  pressure  distributions  at  a  fixed  time,  and  these  data 
are  somewhat  noisy  due  to  the  stairstep  boundary  condition 
modeling  of  the  ramp;  for  RGR  >  30  ft,  peaks  were  evaluated 
by  smoothing  the  pressure  time  histories  two  cells  above 
the  ramp,  giving  a  smooth  pressure-range  curve.*  Note 
that  the  second  peaks  are  in  reasonably  good  agreement  with 


* 

Unfortunately  the  pressure  histories  for  RGR  <  30  ft  were 
not  available  for  data  analysis. 
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Che  prediction  technique  used  to  design  the  simulator.  Also 
note  that  for  RGR  _>  ^0  ft  the  first  and  second  peaks  are 
equal. 

The  calculated  shock  arrival  times  for  the  first  and 
second  static  pressure  peaks  are  included  in  Fig.  10 
The  arrival  time  difference  between  peaks  grows  rapidly  for 
the  first  30  feet  up  the  ramp,  and  then  remains  constant  at 
about  lms/SOOT^^,  In  addition,  Fig.  10  depicts  the  Mach 
stem  growth  versus  ramp  ground  range.  The  top  of  the  Mach 
stem  traces  a  path  at  an  average  angle  of  about  9  degrees 
above  the  ramp  surface,  which  is  consistent  with  shock  tube 
data  for  square  wave  shock  reflections  from  wedges  (Ref.  7). 
Note  that  the  Mach  stem  growth  for  the  equivalent  nuclear 
case  is  more  rapid  than  in  the  case  of  the  simulator. 

Calculated  static  pressure  histories  are  presented 
in  Fig.  11  for  various  stations  on  the  ramp  (34  ft  <  RGR  < 

60  ft).  The  second  peak  dominates  for  RGR  34  ft,  and 
then  gradually  melts  into  backside  of  the  waveform.  For 
RGR  >  60  ft,  the  second  peak  has  essentially  disappeared. 
Comparisons  of  static  pressure  histories  at  h  *  0,  1  and 
5.5  ft  normal  to  the  ramp  for  station  17  indicate  that 
there  is  no  vertical  pressure  gradient  on  the  front  end  of 
the  waveform. 

Fig.  12  gives  the  calculated  dynamic  pressure  his¬ 
tories  on  the  ramp  at  stations  corresponding  to  the  static 
pressure  histories  of  Fig.  11.  At  small  ground  ranges, 
the  second  peak  dominates  the  first  peak.  The  second  peak 
decays  in  magnitude  and  duration  as  the  Mach  stem  pro¬ 
gresses  up  the  ramp,  and  has  essentially  disappeared  for 
RGR  >  60  ft.  Comparisons  of  dynamic  pressure  histories  at 
h  <*  0,  1  and  5.5  ft  normal  to  the  ramp  for  station  17  indi¬ 
cate  very  little  vertical  gradient  for  times  less  than  0.8 
ms  after  shock  arrival.  However,  the  h  *  1  ft  station 
shows  a  strong  second  peak  at  about  1  ms  which  is  absent 
from  the  h  ■  0  and  5  ft  records.  We  believe  that  this  is 
caused  by  a  high  density  slug  of  gas  at  this  altitude.  A 
slip  line  (with  high  density  material  above  and  lower  den¬ 
sity  material  below)  emanates  from  the  triple  point.  As 
the  slip  line  approaches  the  ramp  it  curls  forward  forming 
a  region  of  high  densiry  fluid  near  the  ramp  surface  (h 
1  ft/500T^/^)  while  the  Mach  stem  at  this  station  is  abaut 
10  ft  high.  This  effect  is  similar  to  the  contact  surface 
rollup  observed  in  numerical  simulations  of  nuclear  HOB 
detonations  and  square  wave  shock  reflections  from  wedges 
(Ref.  9).  This  increase  in  dynamic  pressure  near  the  ramp 
can  be  very  important  to  airblast  loads  on  above  ground 


7 


structures — it  increases  both  the  peak  loads  and  the 
impulses  to  approximately  2  ms/SOOT3-/3. 

Lec  us  now  relate  che  simulator  environment  to  an  equi¬ 
valent  nuclear  height-of-burst  case.  Fig.  13  gives  the 
ideal,  nuclear  peak  overpressure  HOB  curves  as  constructed 
by  H.  J.  Carpenter  (Ref.  10).  Region  A  corresponds  to  the 
regular  reflection  regime,  and  region  B  corresponds  to  the 
Mach  reflection  regime  where  the  static  pressure  waveforms 
on  the  ground  contain  two  peaks.  In  regions  and  B2> 
first  and  second  peaks  dominate,  respectively.  Along  the 
dashed  curve  the  first  and  second  peaks  are  equal.  Figure 
9  indicates  that  for  30  ft  <  RGR  <  60  ft,  first  and  second 
peaks  are  equal  and  range  from  600  psi  down  to  400  psi. 
Figure  13  then  indicates  that  for  this  range  in  pressure, 
the  nuclear  HOB  parameters  are  the  following: 

100  ft/KT1/3  _<  HOB  <_  120  ft/KT1,/3 

190  ft/KT1/3  £  OR  <  210  ft/KT1/3 

Thus  the  simulator  as  analyzed  in  this  report  gives  an  air- 
blast  environment  which  is  equivalent  to  a  nuclear  detona¬ 
tion  at  height-of-burst  of  about  110  ft/KT1'3  and  a  ground 
range  of  about  200  ft/KT3  . 

Finally,  let  us  consider  the  effective  yield  of  the 
simulator.  A  500T  high  explosives  surface  burst  produces 
a  blast  wave  flow  field  which  is  equivalent  to  about  a  1-KT 
nuclear  surface  burst  (or  a  2-KT  nuclear  free  air  burst) . 
Nuclear  static  pressure  waveforms  in  the  400  psi  to  600  psi 
Mach  reflection  regime  have  double  peaks  with  a  time  separa¬ 
tion  between  peaks  of  about  2  ms/KTjj3// 3  (Ref.  10).  The 
FAST2D  calculation  of  the  simulator  flow  field  indicates  a 
time  separation  between  peaks  of  about  1  ms/SOOjijgjgg ,  i.e., 
the  time  separation  for  the  simulator  is  too  small  by  a 
factor  of  about  two.  We  believe  that  the  time  separation 
between  peaks  can  be  increased  by  making  the  Mach  stem 
climb  more  rapidly.  This  can  be  accomplished  by  simul¬ 
taneously  decreasing  the  ramp  angle  and  moving  the  ramp 
toward  the  charge. 

V.  SUMMARY  AND  CONCLUSIONS 

Height-of-burst  detonations  create  airblast  environments 
and  diffraction  loads  which  are  more  severe  than  the  surface 
burst  c.  se  in  the  high  overpressure  Mach  reflection  regime. 
There  is  an  ongoing  need  to  simulate  these  HOB  environments 
on  a  large  scale  to  validate  the  survivability  of  military 
systems  to  blast  effects.  We  propose  using  an  existing 
high  explosives  test  bed,  say  a  SOOT  hemispherical  charge. 
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co  create  the  iree  field  blast  environment.  A  large  ramp 
would  be  located  near  the  charge.  Shock  diffraction  on  the 
ramp  generates,  in  a  rather  natural  wav,  a  flow  field  which 
simulates  the  HOB  blast  environment  with  high  fidelity. 

A  parametric  analysis  of  such  HOB  simulators  indicates 
that  a  30°  ramp  situated  about  200  feet  from  a  SOOT  hemis¬ 
pherical  charge  would  give  useful  environments.  The  flow 
field  details  near  such  a  ramp  were  investigated  with  a  2-D 
hydrocode  calculation.  The  calculation  indicates  that 
double  peaked  static  and  dynamic  pressure  waveforms  were 
created  near  the  ramp  surface.  In  the  400  to  600  psi  range, 
the  calculated  first  and  second  static  pressure  peaks  were 
equal.  By  use  of  the  nuclear  HOB  curves,  it  was  determined 
that  the  blast  flow  field  corresponds  to  a  nuclear  detona¬ 
tion  at  a  height-of-burst  of  100  to  120  ft/KTN1''3  and  a 
ground  range  of  190  to  210  ft/KT{jl/3.  Time  separations 
between  static  pressure  peaks  were  found  to  be  about  1  ms/ 
500T^^gg;  this  value  was  too  small  by  about  a  factor  of 
two  for  the  nuclear  case. 

VI.  RECOMMENDATIONS 

Additional  analysis  should  be  performed  to  refine  the 
HOB  simulator  design.  The  2-D  hydrocode  simulations  are 
quite  useful  because  they  allow  one  to  examine  the  entire 
flow  field  in  a  non- interfering  way.  An  improvement  is 
needed  on  the  boundary  condition  modeling  of  the  ramp — the 
stairstep  model  gave  very  noisy  results  on  the  ramp  surface. 
Small  charge  (say  4-lb  hemispherical  PBX-9404  charges)  tests 
can  provide  an  experimental  definition  of  the  blast  environ¬ 
ment.  Ramp  angle,  location  and  surface  curvature  could  be 
varied  parametrically  in  such  tests.  Pressure  gauges  on 
the  ramps  can  measure  static  pressure  histories  with  high 
fidelity,  while  shadowgraph  photography  can  capture  the 
shock  structure  on  the  ramp.  These  results  could  be  used 
to  design  a  simulator  which,  we  suggest,  should  be  fielded 
on  the  next  500T  HE  test. 
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Fig.  2  —  HOB  simulator  concept  (graded  ramp)  on  a  large-scale  HE  test 
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Fig.  3  —  Comparison  of  calculated  pressure  and  density  contours  for  a  nuclear 
HOB  case  and  a  square  wave  shock  on  a  wedge 
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Fig.  5  —  Parametric  results  of  peak  reflected  pressures  on 
the  ramp  HOB  simulator 


Oiapman-Jouguet  initial  conditions  Results 


Fig.  7  —  Calculated  pressure  and  density  contours  at  times  At  =  3.0,  5.61, 
9.27,  and  13.4  ms/500  T1'3  (t0  =  19.2  ms/500  T1'3) 
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Fig.  8  —  Calculated  flow  details  at  At  =  9.61  ms/500  T1/® 
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Fig.  v>  —  Shadowgraph  of  the  shock  wave  structure  formed  by  an  8-lb  TNT-driven 
blast  wave  i  HOB  =  1.04  ft/lb1'®)  diffracting  on  a  31s  ramp:  incident  pressure  at 
the  beginning  of  the  ramp  was  about  120  psi.  (Courtesy  of  W.  F.  Duaziak. 
Information  Science.  Inc.) 
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Fig.  12  —  Calculated  dynamic  t  <tssurt 
the  ramp  (tc  =  19 
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Fig.  13  —  Ideal  nuclear  poak  overpresaure  height-of-burat  curves  (Ref.  10) 
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